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We present high precision estimates of the exponents of a quan- 
tum phase transition in a planar antiferromagnet. This has been 
made possible by the recent development of cluster algorithms for 
quantum spin systems, the loop algorithms. Our results support 
the conjecture that the quantum Heisenberg antiferromagnet is in 
the same universality class as the 0(3) nonlinear sigma model. The 
Berry phase in the Heisenbrg antiferromagnet do not seem to be 
relevant for the critical behavior. 

I. INTRODUCTION 

Instead of classical transitions controlled by temperature T a quantum phase 
transition between a symmetry broken phase with long-range Neel order and 
a quantum disordered state with a finite spin excitation gap may be realized 
at T = by controlling a parameter g to increase quantum fluctuations. Crit- 
icalities around such quantum phase transitions at g = g c may reflect inherent 
quantum dynamics of the system and yield unusual universality classes with 
rich physical phenomena. 

The most prominent example are the high temperature superconductors. 
There the quantum spin fluctuations are thought to lead to d-wave supercon- 
ductivity as soon as antifcrromagnetism is suppressed by hole doping. This 
close connection between antiferromagnetism, quantum fluctuations and high 
temperature superconductivity has triggered many theoretical investigations. 

Most of these investigations are based on a mapping to an effective field 
theory, the 2D 0(3) quantum nonlinear sigma model. This sigma model is 
in the same universality class as the 3D 0(3) classical sigma model or the 
3D classical Heisenberg model. A large number of detailed predictions about 
quantum critical behavior has been made for the sigma model [jl],^j. However 
the spin- 1/2 quantum antiferromagnet generally contains Berry phase terms 
H that are not present in the sigma model. The relevancy of these terms is 
not clear. 

In order to shed light onto this question we have simulated a two dimen- 
sional quantum antiferromagnet (2D QAFM) that exhibits a quantum phase 
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transitions. We calculate the critical exponents to determine the universality 
class and to check predictions made based on the nonlinear sigma model. First 
results have been published in Ref. (J]. 



A. Quantum critical exponents 

The critical exponents of a quantum phase transition at T — can be de- 
fined similar to a classical finite temperature phase transition. The quantum 
mechanical control parameter g plays the role of the temperature in the clas- 
sical system. Approaching the quantum critical point from the disordered side 
(<? > g c ) the correlation length diverges as 

Z<x(g-g c )-". (l) 

By the Trotter-Suzuki mapping the <i-dimensional quantum system can be 
mapped onto a d + 1-dimensional classical system. At zero temperature the 
system is infinite also in the additional imaginary time direction. The space 
and time dimensions are however not necessarily equivalent, and the correlation 
length in the time direction diverges in general with a different exponent 

eraCff-ffc)-*", (2) 

where z is the dynamical exponent. In a Lorentz invariant system space and 
time directions are equivalent and z = 1. Related to the divergence of the 
correlation length is a vanishing of the spin excitation gap 

A oc (g - gc ) z r (3) 

When passing through the critical point long range order is established. The 
order parameter in the case of a Neel ordered antiferromagnet is the staggered 
magnetization 
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(4) 



where N is the number of spins in the lattice, S* the z-component of the spin at 
site r and Q = (tt,tt). Close to the critical point the staggered magnetization 
behaves as 

m s oc (g c - g) p . (5) 
At the critical point itself the real space correlation show a power-law falloff 

(5^ r z )oce ,Qr r- (d+2 - 2+,,) , (6) 

where rj is the correlation exponent. These three exponents are related by the 
usual scaling law 

2f3=(d+z-2 + r))v, (7) 
where the effective dimension is d + z in a quantum system. 
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TABLE I. Critical exponents /3, v, and r\. Listed are both the estimates without 
making any assumption for z, and the best estimate if Lorentz invariance (z — 1) is 
assumed. For comparison the exponents of the 3D classical Heisenberg (0(3)) model, 
the 3D Ising model and the 2D quantum mean field exponents are listed. 



model 


V 





V 


2D QAFM 

Lorentz invariant 2D QAFM 


0.685 ± 0.035 
0.695 ± 0.030 


0.345 ± 0.025 
0.345 ± 0.025 


0.015 ±0.020 
0.033 ± 0.005 


3D 0(3) I 
3D Ising | 
mean field 




0.7048 ± 0.0030 
0.6294 ± 0.0002 
1 


0.3639 ± 0.0035 
0.326 ± 0.004 
1/2 


0.034 ± 0.005 
0.0327 ± 0.003 




TABLE II. Universal prefactor fh(oo) in the linear temperature dependence of the 
uniform susceptibility at criticality. Listed are the results for the quantum nonlinear 
sigma model in a 1/7V expansion, the results by classical Monte Carlo simulation on 
a 3D classical rotor model and the result of the present study. 



method 


Ref. 


Qi (oo) 


1/N expansion 
classical Monte Carlo 

quantum Monte Carlo t 


2 
2 
b 


is study 


0.2718 
0.25 ± 0.04 
0.26 ±0.01 



B. Predictions from the nonlinear sigma model 

As mentioned above most analytic calculations of quantum critical behavior 
are based on the 2D 0(3) quantum nonlinear sigma model (QNLcrM). Here 
we want to review the critical properties of the sigma model relevant for the 
current study. 

The critical exponents of the QNLcrM can be determined from simple sym- 
metry, universality and scaling arguments jl]]^] . Lorentz invariance implies that 
z = 1. Furthermore the 2D QNLcrM is equivalent to the 3D classical sigma 
model. This in turn is in the universality class of the 3D classical 0(3) model, 
or the classical 3D Heisenberg ferromagnet. The exponents /?, v and r\ should 
thus be the same as the well known classical exponents of these models (see 
Tab. I). 

Chakravarty, Halperin and Nelson have discussed the phase diagram of a 
planar Heisenberg antiferromagnet in the framework of the QNLcrM. They con- 
centrate on the ordered phase and describe it as a classical 2D antiferromagnet 
with renormalized parameters. 

Chubukov, Sachdev and Ye have investigated the quantum critical regime of 
the QNLcrM in close detail They make some further predictions based on 
scaling arguments. On the ordered side the spin stiffness p s vanishes as 

Ps K(g c -gY d +*-^ = (g c - g y, (8) 
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where the second equivalence conies from the prediction that 2 = 1. Addition- 
ally it follows from general scaling arguments that the uniform susceptibility 
at the critical point is universal: 




Here c is the spin wave velocity and fii(oo) a universal constant. Estimates for 
fii(oo) are listed in Tab. II. 

The spin wave velocity c scales as 

(10) 

and is thus regular at the critical point if z = 1. 



C. What about Berry phases? 

The equivalence of the 2D QAFM to the 2D QNLcrM however is still an open 
question because of the existence of Berry phase terms in the QAFM that are 
not present in the QNLcrM || . It has been argued that these terms cancel in 
special cases, such as in the bilayer model [RB. Then it is plausible that the 
quantum phase transition is in the same universality class as the QNLcrM. This 
was confirmed by quantum Monte Carlo calculations of Sandvik and coworkers 
J7|,p|jl0|] ■ They have investigated the finite size scaling of the ground state 
structure factor and susceptibilities on lattices with up to 10 x 10 x 2 spins. 
Although these lattices are quite small they still found good agreement of the 
exponents z and r\ with the QNLcrM predictions |^,^| . 

In another study Sandvik et al. JlC| ] have investigated finite temperature 
properties of the bilayer QAFM on larger lattices and also found good agree- 
ment with the QNLcrM predictions. In the absence of Berry phase terms the 
equivalence of the QAFM and the QNLcrM is quite well established by these 
simulations. 

But in general these Berry phase terms exist. Chakravarty et al. argue that 
they can change the critical behavior and lead to different exponents 
Chubukov et al. on the other hand argue that the Berry phase terms are 
dangerously irrelevant Q and do not influence the critical behavior. 

Previous numerical simulations on dimerized square lattices [p|Jl2| are indeed 
not consistent with the QNLcrM predictions. Sandvik and Vekic || find a 
dynamical exponent z ^ 1, but their largest system was only 10 x 10 spins. 
The deviation could be a problem with scaling arising from inequivalent spatial 
directions. 

Katoh and Imada however found j»l, compatible with Lorentz invariance. 
Additionally they calculated the correlation length £ and from it the exponent 
v. The validity of their result v w 1 is however again questionable because 
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because of the restriction to very small lattices of 12 x 12. On the other hand 
the discrepancy could be an effect of the Berry phase terms that are present in 
the dimerized square lattice but probably not in the bilayer. 

The main purpose of the simulations reported is to she light onto this question 
and to clarify the role of the berry phase terms. Our results support the ideas 
of Chubukov et. al. Q that the Berry phase terms are dangerously irrelevant. 

II. ALGORITHM AND PARALLELIZATION 

Using the new quantum cluster algorithms, the loop algorithms [ p^|Jl4| it 
is possible to simulate much larger lattices at lower temperatures, just as the 
corresponding classical cluster algorithms have allowed the simulation of critical 
classical spin systems. With these algorithms it has for the first time become 
possible to study quantum critical spin systems in detail. 

A disadvantage of the cluster methods however is that they cannot be vector- 
ized as easily as the local update algorithms. Using powerful vector machines 
is therefore not an option. Fortunately however most of the modern super- 
computers are parallel machines, and Monte Carlo methods are nearly ideally 
suited for that architecture. 

One of the authors has developed an object oriented Monte Carlo library 
in C++ Using this library it is very simple to parallelize a Monte Carlo 
program and to port it to new parallel computers. The library automatically 
parallelizes any Monte Carlo simulation at the two "embarrassingly parallel" 
levels. The first level of trivial parallelism is the parameter parallelism. Simu- 
lation with different parameters, such as system size, coupling or temperature 
can be performed independently in parallel. At this level there is practically 
no overhead due to the parallclization. We get perfect speedup and the library 
takes care of load balancing. 

A single simulation can similarly be parallelized by running it in parallel 
with different initial states and random seeds on each of the processors. The 
simulations run nearly independent. Communication is required only at the 
start and the end of the simulation. This level of parallelization incurs some 
overhead however. The overhead is the time used to thermalize a simulation. 
We loose efficiency if this thermalization time becomes comparable to the time 
actually needed for the simulation. 

The third and deepest level of parallelization cannot be automatically done 
by the library since it depends on the algorithm used for Monte Carlo. The 
lattice used for one simulation can be spread over many processors. This par- 
allclization has to be done by the programmer of the algorithm, but it is sup- 
ported by various functions of the library. It is worthwhile to invest time in 
this parallelization only in two cases. The first is when, as mentioned above, 
thermalization is slow. Often the main reason is however different one. Large 
lattices simply might not fit into the memory of one processor. 
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FIG. 1. Lattice structure of the 1/5-th depleted square lattice of CaV^g. The 
dashed square indicates the eight spin unit cell used in our calculations. 

In our simulations reported here we have used the 1024-node, 300 GFlop 
Hitachi SR2201 massively parallel computer of the university of Tokyo. At the 
time of its introduction this machine was the fastest general purpose computer 
in the world. Each processor has 256 MByte of local memory, enough to simu- 
late quantum spin systems with 20000 spins at temperatures as low as T = 0.01. 
This was large enough for the present study and we did not spend time on the 
third level of parallclization but used only the first two levels provided by the 
library. 

The algorithm used was the continuous time loop algorithm flifl . The loop 
algorithms, first developed by Evertz et al. |jl3f are quantum version of the 
classical cluster algorithms. The continuous time version is preferable over the 
earlier discrete time versions since it eliminates the need to extrapolate in the 
finite Trotter time step At. In our experience we found that this leads to 
a four-fold speed increase. Additionally the continuous time algorithm uses 
only 10% of the memory compared to the discrete time algorithm, allowing the 
simulation of larger lattices. 

III. THE CAV4O9 LATTICE 

As the universality class of a phase transition does not depend on the micro- 
scopic details of the lattice structure we are free to choose the best lattice for 
our purposes. We have chosen the CaV40g lattice, a 1/5-th depleted square 
lattice depicted in Fig. [I] for our calculations. There are three reasons for this 
choice. Firstly the Berry phase terms are present on this lattice jl|]. Next 



G 




FIG. 2. Phase dia gram of the CaV40g spin lattice as a function of the ratio Jo/ Ji, 
reprinted from Ref. p7| . The leftmost point corresponds to the dimer limit Jo = 
and the rightmost point to the plaquette limit Ji = 0. Circles indicate quantum 
Monte Carlo results for the spin gap, normalized by Jo + Ji. Diamonds show the 
staggered magnetization. As reference the perturbation theory estimates for the gap 
| p^| and the linear spin wave theory (LSW) estimates for the staggered moment have 
been included. 



both space directions are equivalent, in contrast to the dimerized square lattice 
P,^2|. This makes the scaling analysis easier. Finally at the quantum critical 
point all the couplings are nearly equal in magnitude, which is also optimal 
from a numerical point of view. We have performed our simulations on square 
lattices with N = 8n 2 spins, where n is an integer. Our largest lattices con- 
tained 20 000 spins. For the following discussion it is useful to introduce the 
linear system size L in units of the bond lengths a of the original square lattice: 
L = A/5iV/4a. 

The phase diagram of this lattice has been discussed in detail in Ref. jL7j and 
is shown in Fig. @. By removing every fifth spin we obtain a lattice consisting 
of four-spin plaquettes linked by dimer bonds. We label the couplings in a 
plaquette Jo and the inter-plaquette couplings J\. By controlling the ratio of 
these couplings J\/Jq we can tune from Neel order at Ji = Jo to a quantum 
disordered "plaquette RVB" ground state with a spin gap A = J at J\ = 0. 
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FIG. 3. Plot of the ratio of correlation length divided by system size £l/L. At the 
critical point the correlation length calculated in a finite system is proportional to 
the system size. This is the case for (Ji/Jo) c = 0.939 ± 0.001. 

IV. RESULTS 

A. The critical point 

The first step in the determination of the critical behavior is a high precision 
estimate of the critical coupling ratio Ji/Jq. We have calculated the second 
moment correlation length £l on systems of various sizes L. This can be de- 
termined in the usual way from the magnetic structure factor S(q) close to the 
Neel peak at Q: 

g(Q+ ^ = i + awio a +0( ^ q4) - 

The temperature was chosen to be T — Jq/L, keeping the finite 2 + 1 dimen- 
sional system in the cubic regime. From standard finite size scaling arguments 
it follows that this correlation length £l scales proportional to the system size 
L at criticality. We have calculated the ratio £l/L (shown in Fig. ||) for a 
variety of couplings and system sizes up to N = 9600. Independence of the 
system size was seen at the critical coupling ratio {Jx/Jq) c — 0.939 ± 0.001. 
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FIG. 4. Finite size scaling of the staggered structure factor and susceptibility at 
the critical point. 



B. The exponents 

Next we have calculated the finite size scaling of both the staggered structure 
factor 5(Q) = L 2 m s and of the corresponding staggered susceptibility. At 
criticality they scale like 

5(Q) oc L 2 - Z ~ r > (12) 
Xs oc L 2 - r ' (13) 

The temperature was chosen to be T — Jq/(6L), low enough to see ground 
state properties within our accuracy. By fitting the results shown in Fig. ^ we 
obtain the estimates z = 1.018 ± 0.02 and r\ = 0.015 ± 0.020. This is perfectly 
consistent with the Lorentz invariance (z — 1) expected from a mapping to the 
QNLcrM. We will discuss r\ below together with the other exponents. From 
these fits it is also obvious that at least N = 800 spins are necessary to obtain 
good scaling. 

The remaining exponents [3 and v are best calculated from the magnetization 
m s and the spin stiffness p s on the ordered side. Good estimates for m s and p s 
can be obtained by the Hasenfratz-Niedermayer equations jl9| . These authors 
have calculated the exact finite-size and finite-temperature values of the low- 
temperature uniform and staggered susceptibilities \ u and Xs for the ordered 
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FIG. 5. Staggered magnetization m s and spin stiffness p s calculated by a fit of 
the low temperature susceptibilities on finite lattices to the Hasenfratz-Niedermayer 
equations. 



phase of a 2D QAFM on a lattice with the symmetries of a square lattice. Their 
equations for the staggered susceptibility 



Xs(T,L) 



2r2 



MjL 
3T 



1 



p s Ll 



he 

p s Ll 



[/?i(0 2 + 3/3 2 (0] 



(14) 



and for the uniform susceptibility 
2 Ps 



Xu(T,L) 



1 i,c m- 



3(hc) 2 { 3 p s Ll 



1 



Tie 



3 \p s Ll 



(15) 



are correct for the low temperature regime fcgT <C 27rp s with cubic geometry 
I 3 = « 1. Up to second order in T (or 1/L respectively) the susceptibilities 
are universal, determined by only three parameters: the staggered magnetiza- 
tion m s , the spin stiffness p s and the spin wave velocity c. The shape functions 
/32) /3i; @2 and ip are known exactly for square lattice geometries. Two high 
precision quantum Monte Carlo studies have confirmed the validity of these 
equations for the square lattice QAFM |14] , po[ . 
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We have calculated the susceptibilities for a wide range of couplings 0.95 < 
Ji/Jo < 1-1, lattice sizes 800 < N < 16200 and temperatures 0.006 < 
T/Jo < 0.1. The fits to the Hasenfratz-Niedermayer equations are all ex- 
cellent, with % 2 /d.o.f. ~ 1.5. This is another confirmation of the universality 
of the Hasenfratz-Niedermayer equations. From the fits we obtain the stag- 
gered magnetization m s , the spin stiffness p s and the spin wave velocity c. The 
exponents (3 and v can then be obtained in a straightforward way (see Fig. |^) 
and are listed in Tab. I. 

C. Discussion 

Let us now discuss the results. First we observe that the exponents satisfy 
the scaling relation Eq. (0) , indicating the validity of the scaling ansatz for this 
quantum phase transition. The exponents (3, v and r\ are in excellent agreement 
with the exponents of the 3D classical 0(3) or Heisenberg model. They are 
however incompatible with the mean field exponents calculated by Katoh and 
Imada on small lattices. 

Assuming Lorentz invariance z = 1 we can improve our estimates for the 
other three exponents. The agreement of the improved estimates with the 
3D 0(3) exponents becomes even better. We can also rule out the 3D Ising 
universality class whose exponents are also listed in Table I for a comparison. 

This excellent agreement is a strong indication that the Berry phase terms 
in the 2D QAFM are indeed dangerously irrelevant as suggested by Chubukov, 
Sachdev and Ye. To further confirm their predictions we have calculated the 
uniform susceptibility close to criticality down to T = 0.02, more than an order 
of magnitude lower than Ref. m . We have extrapolated the finite size results on 
lattices with up to N = 20000 spins to the thermodynamic limit. Looking for 
the coupling at which a linear behavior occurs gives an independent estimate 
of the critical point: (Ji/Jo) c = 0.939 ± 0.002, in excellent agreement with 
the above estimate. The linear slope is Hi (oo) (a Jo /fi-c) 2 = 0.238 ± 0.003. By 
extrapolating the spin wave velocity determined in the ordered phase by the 
Hasenfratz-Niedermayer fit to the critical point we get Hc/aJo = 1.04 ± 0.02 
and thus f2i(oo) = 0.26 ±0.01, again in excellent agreement with Chubukov et 
al. (see Tab. II). 

V. SUMMARY AND OUTLOOK 

To summarize, we have performed a large scale quantum Monte Carlo sim- 
ulation of a quantum phase transition in a planar antiferromagnet. The new 
quantum cluster algorithms, in particular the continuous time loop algorithm 
allow high precision simulation of critical quantum systems. 
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The critical exponents that we have calculated (listed in Table I) agree within 
our errors with the exponents of the classical 3D 0(3) or Heisenberg model. 
The dynamical exponent z = 1.018 ± 0.02, consistent with Lorentz invariance. 
This is compelling numerical evidence for the conjecture that the quantum 
Heisenberg antiferromagnet is in the same universality class as the 2D quantum 
nonlinear sigma model and the 3D Heisenberg ferromagnet. The Berry phase 
terms that are present in the Heisenberg antiferromagnet for non-integer spin 
do not seem to be influence the critical behavior. This supports the conjecture 
by Chubukov, Sachdev and Ye || that they are dangerously irrelevant. 

While the accuracy achieved in the present simulation is remarkable for a 
simulation of a quantum system it is not very good compared to the best 
classical results that are an order of magnitude more accurate. If one wishes for 
higher accuracy the best approach could be to generalize the histogram methods 
to quantum systems and to do a finite size scaling study of cumulants, similar to 
the ones done in the classical case || . Even then we might not be able to reach 
the same accuracy as in the classical simulations for two reasons. Firstly we 
cannot measure the order parameter m s in a quantum Monte Carlo simulation, 
but only its square m 2 s . Therefore we can only calculate the cumulants of 
rrig, which will be less favorable numerically. Another difference between a 
classical and a quantum system is that in the 3D classical system all three 
space directions are equivalent. In the (2+l)D quantum system on the other 
hand the time direction is not equivalent to the space direction, which makes 
a scaling analysis more complex. 
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